Spinal fMRI demonstrates segmental organisation of functionally connected networks in the cervical spinal cord: A test–retest reliability study

Abstract Resting functional magnetic resonance imaging (fMRI) studies have identified intrinsic spinal cord activity, which forms organised motor (ventral) and sensory (dorsal) resting‐state networks. However, to facilitate the use of spinal fMRI in, for example, clinical studies, it is crucial to first assess the reliability of the method, particularly given the unique anatomical, physiological, and methodological challenges associated with acquiring the data. Here, we characterise functional connectivity relationships in the cervical cord and assess their between‐session test–retest reliability in 23 young healthy volunteers. Resting‐state networks were estimated in two ways (1) by estimating seed‐to‐voxel connectivity maps and (2) by calculating seed‐to‐seed correlations. Seed regions corresponded to the four grey matter horns (ventral/dorsal and left/right) of C5–C8 segmental levels. Test–retest reliability was assessed using the intraclass correlation coefficient. Spatial overlap of clusters derived from seed‐to‐voxel analysis between sessions was examined using Dice coefficients. Following seed‐to‐voxel analysis, we observed distinct unilateral dorsal and ventral organisation of cervical spinal resting‐state networks that was largely confined in the rostro–caudal extent to each spinal segmental level, with more sparse connections observed between segments. Additionally, strongest correlations were observed between within‐segment ipsilateral dorsal–ventral connections, followed by within‐segment dorso–dorsal and ventro–ventral connections. Test–retest reliability of these networks was mixed. Reliability was poor when assessed on a voxelwise level, with more promising indications of reliability when examining the average signal within clusters. Reliability of correlation strength between seeds was highly variable, with the highest reliability achieved in ipsilateral dorsal–ventral and dorso‐dorsal/ventro–ventral connectivity. However, the spatial overlap of networks between sessions was excellent. We demonstrate that while test–retest reliability of cervical spinal resting‐state networks is mixed, their spatial extent is similar across sessions, suggesting that these networks are characterised by a consistent spatial representation over time.

reliability achieved in ipsilateral dorsal-ventral and dorso-dorsal/ventro-ventral connectivity.However, the spatial overlap of networks between sessions was excellent.
We demonstrate that while test-retest reliability of cervical spinal resting-state networks is mixed, their spatial extent is similar across sessions, suggesting that these networks are characterised by a consistent spatial representation over time.

Practitioner Points
• We studied spinal resting-state networks and their test-retest reliability.
• Dorsal and ventral networks emerged within segmental levels and strong within-segment connections were seen across grey matter horns.
• Reliability estimates were mixed but the spatial overlap of networks was excellent.

| INTRODUCTION
Spinal cord functional magnetic resonance imaging (fMRI) is a novel but rapidly developing field (Kinany, Pirondini, Micera, & Van De Ville, 2022;Powers et al., 2018).Combined with brain fMRI, it holds promise for investigation of information processing across all levels of the central nervous system in both health and disease.
Like the brain, the spinal cord is characterised by spontaneous fluctuations in the blood oxygen level dependent (BOLD) signal in the absence of overt stimulation.This intrinsic activity of the spinal cord has been shown to form organised resting-state networks, which can be broadly divided into motor and sensory (Harrison et al., 2021).
Reports of strong temporal correlations between the sensory (dorsal) horns and motor (ventral) horns within the cervical spinal cord have dominated the spinal fMRI resting-state literature (Barry et al., 2014(Barry et al., , 2016;;Eippert et al., 2017;San Emeterio Nateras et al., 2016;Weber et al., 2018).Furthermore, unilateral sensory networks have also been observed in resting spinal data, which were limited in rostro-caudal extent, corresponding to the underlying segmental anatomy of the cord (Kong et al., 2014).Early evidence from simultaneous brain-spine fMRI has also shown that spinal and cerebral resting-state networks are correlated, suggesting a unified functional architecture of intrinsic networks in the central nervous system (Vahdat et al., 2020).
Brain resting-state fMRI is frequently used as a biomarker for identification of neurodivergent states/conditions or treatment effects (Drysdale et al., 2017;Pfannmöller & Lotze, 2019;Taylor et al., 2021).Reliable detection of resting-state networks in the spine would extend this approach to information processing occurring at the level of the cord, such as early modulation of noxious signals or motor functioning (Kinany, Pirondini, Micera, & Van De Ville, 2022;Tinnermann et al., 2021).Acquiring fMRI recordings from the spinal cord, however, faces unique anatomical, physiological, and methodological challenges, including, among others, the small size of the cord, influence of physiological noise, and reliable static magnetic field shimming (Kinany, Pirondini, Micera, & Van De Ville, 2022;Tinnermann et al., 2021).These challenges can limit the quality of obtained data and thus pose a threat to the reliability of spinal fMRI.To date, the few studies that investigated the reliability of resting-state spinal cord fMRI showed good test-retest reliability (intraclass correlation coefficient [ICC] = 0.64-0.7) in network properties using graph theory measures at 3 T (Liu et al., 2016) and fair reliability (ICC = 0.54-0.56) in region-to-region connections at 7 T (Barry et al., 2016).A recent assessment of reliability of region-to-region connections at 3 T has further shown that reliability was fair to good for dorso-dorsal and ventro-ventral connections but poor for within and betweenhemicord connections across the cervical cord and generally poor for all connections within individual segmental levels (Kaptan et al., 2022).These studies, however, assessed test-retest reliability within the same scanning session.Given that longer lag between scans is associated with poorer reliability in cerebral fMRI (Bennett & Miller, 2010, 2013) and that the scanning setup for spinal cord fMRI is considerably more complicated than for cerebral fMRI (Kinany, Pirondini, Micera, & Van De Ville, 2022;Powers et al., 2018;Tinnermann et al., 2021), investigations of test-retest reliability of spinal cord fMRI that span separate scanning sessions are warranted.Such investigations will indicate the feasibility of using spinal cord fMRI to reliably detect the effects of experimental manipulation or clinical interventions across different visits, such as perturbations related to experimental pain, persistent pain (e.g., post-surgical), or treatment effects.
Test-retest reliability is inherently tied to data quality.Acquiring good quality spinal cord fMRI recordings is complicated by the influences of baseline physiology and susceptibility artefacts related to differing magnetic susceptibility profiles of surrounding tissues (Kinany, Pirondini, Micera, & Van De Ville, 2022;Saritas et al., 2014;Tinnermann et al., 2021).Shimming procedures can minimise the effects of these factors by reducing magnetic field inhomogeneities.A combination of high order and slice-specific z-shimming is frequently used in spinal cord fMRI to improve signal quality (Eippert et al., 2017;Finsterbusch et al., 2012;Kinany, Pirondini, Mattera, et al., 2022;Vahdat et al., 2020).Nonetheless, while z-shimming offers large signal gains by accounting for the off-resonance variation along the cord, implementing simultaneous x, y, and z-shimming can achieve additional benefits by preventing signal loss caused by magnetic field gradients in left/right and anterior/posterior directions (Alonso-Ortiz et al., 2023;Islam et al., 2019).Furthermore, given that magnetic field inhomogeneities can induce artefacts in traditional echo planar imaging (EPI) sequences incorporating fat saturation pulses, using a spectral-spatial pulse exciting only tissue water could further improve signal quality (Bernstein et al., 2004), along with allowing for shorter repetition time (TR) or larger number of slices acquired within the same time.
This study assesses the test-retest reliability of cervical spinal cord resting-state fMRI over two separate scanning sessions.Additionally, we demonstrate a novel implementation for acquiring BOLDsensitive resting-state spinal fMRI and characterise functional connectivity relationships in the cervical cord in healthy adult volunteers.In particular, the acquisition sequence used here operates on a General Electric (GE) scanner platform, using second-order shimming and x, y, and z slice-specific linear shimming, together with spectral-spatial excitation pulses designed to excite tissue water only.This approach reduces signal dropout and increases temporal signal-to-noise ratio (tSNR) within the cervical spinal cord (see Tsivaka et al. (2023) for full details of the acquisition method).
Our pre-registered hypotheses (Kowalczyk et al., 2021) are: 1. Discrete resting-state sensory and motor networks should be observable in regions of the dorsal and ventral cervical spinal cord, respectively, using T2*-weighted BOLD EPI.
2. Spinal responses observed during the assessments of hypothesis 1 will be reliable, with ICC inter-session test-retest reliability statistics greater than 0.4.

| Participants
Data from 23 healthy right-handed (as assessed by the Edinburgh Handedness Inventory; Oldfield, 1971) adult volunteers (13 females, mean + SD age: 23.91 ± 3.84 years) were collected for all study visits and survived all quality assessments.Full details of participant/data exclusion are shown in Figure 1.
Full inclusion and exclusion criteria for this study are outlined in the study preregistration (Kowalczyk et al., 2021).Briefly, participants were excluded due to: (1) history of psychiatric, medical, or psychological conditions, (2) history of substance or alcohol abuse, (3) regular use of medications affecting the central nervous system, (4) irregular menstrual cycle for female participants, (5) MRI-related contraindications.Additionally, participants were excluded if they were unwilling to adhere to the following lifestyle guidelines before each visit: (1) abstain from alcohol for 24 h, (2) limit caffeine consumption to one caffeinated drink on each study day, (3) abstain from non-steroidal anti-inflammatory drugs (NSAIDs) or paracetamol for 12 h, (4) abstain from nicotine-containing products for 4 h.
Additional exclusion criteria based on data quality were used.
Subjects with incomplete functional acquisition, poor image contrast, and low tSNR were excluded (see Section 2.5 for details).
Written informed consent was obtained.This study was approved by the Psychiatry, Nursing, and Midwifery Research Ethics subcommittee at King's College London, UK (HR-16/17-4769).

| Procedure
This study comprised three visits-a screening/familiarisation visit and two identical MRI visits for test-retest purposes.The mean (±SD, range) interval between each study visit was 21 (±22, 1-84) days; mean (±SD, range) interval between scanning visits was 31 (±24, 1-84).Additional assessments not described here, namely evoked response pain and motor fMRI, were collected during the study visits and are detailed in the preregistration (Kowalczyk et al., 2021).

| Session 0: Screening and familiarisation
Compliance with study lifestyle guidelines (see Section 2.1) was assessed at the beginning of the session.Participants underwent breath alcohol and urine drugs of abuse tests to check alcohol/substance use.Caffeine, nicotine, and NSAIDs/paracetamol intake were assessed by self-report.Participants were familiarised with the scanner environment by visiting a mock scanner.

| Sessions 1 and 2: MRI scanning
Sessions 1 and 2 were identical.The sessions began with an assessment of compliance with the study lifestyle guidelines as described above.Additionally, participants completed the state version of the State Trait Anxiety Inventory (STAI; Spielberger et al., 1971) to assess differences in anxiety levels between sessions.No differences were observed (t[22] = 1.23, p = .223,d = 6.12, 95% CI [À0.67; 1.6]; Session 1 mean ± SD = 27.61 ± 1.32; session 2 mean ± SD = 29.17± 1.49).Subsequently, MRI data was collected in the following order: (1) optimisation of static 0th, 1st, and 2nd order shims and linear slicespecific shims, (2) structural data acquisition, and (3) 10 min 50 s resting-state scan (see Section 2.3).Participants were instructed to keep their eyes open and look at the fixation cross displayed in the centre of the screen (white cross on a black background).Respiratory and cardiac traces were recorded with respiratory bellows and a pulse oximeter, respectively, along with scanner triggers (at the start of each TR), throughout the scan.

| MRI acquisition
Data were acquired using a 3 T GE MR750 System (General Electric, Chicago, Illinois) equipped with both a 12-channel head, neck, and spine coil and a 4-channel neurovascular array at the NIHR Wellcome King's Clinical Research Facility, King's College London.A sagittal 3D CUBE T2-weighted structural image was acquired at the beginning of the scanning session over 64 slices with a coverage of the whole brain Four dummy scans were acquired to enable the signal to reach steady-state, followed by 256 volumes.Full details of the acquisition sequence can be found in Tsivaka et al. (2023).For 13 participants, F I G U R E 1 Selection of participants fulfilling the study eligibility criteria and data quality assurance.MR, magnetic resonance, tSNR, temporal signal-to-noise ratio.the manufacturer's EPI internal reference option was used.The internal reference acquires four non-phase-encoded echoes before the EPI echo train, which are used to apply a phase correction to the EPI data.
Upon further inspection of the data, this was shown to contribute to slice misalignment in 5 of the 13 participants (anterior-posterior direction) and thus the setting was disabled for subsequently recruited participants.In order to keep the two MRI visits identical, however, the internal reference was used on both MRI visits for the first 13 participants even after the issue was discovered.
A group mean EPI image in Polytechnique Aix-Marseille University and Montreal Neurological Institute 50 (PAM50) template space (De Leener et al., 2018) showing the average signal intensity (Figure S2), as well as a representative single-subject temporal mean EPI in native space are available for download from NeuroVault (https://identifiers.org/neurovault.collection:13616).

| Data pre-processing
Data were processed using Spinal Cord Toolbox (SCT) version 5.4 (De Leener et al., 2017), AFNI's 3dWarpDrive (Cox, 1996;Cox &Hyde, 1997), andFSL version 6.0.4 (Jenkinson et al., 2012;Smith et al., 2004).Visual quality assurance was performed on raw data and at each stage of processing.Five scans acquired with an early version of the functional sequence using the internal reference (see above) had several slices come out of alignment with the rest of the spinal cord due to a shift in the anterior-posterior (EPI phase-encoding) axis.
A custom in-house Matlab version 9.5.0 (Mathworks Inc.) script was used to move the slices back into alignment with the rest of the cord.
Briefly, for each slice, a 1D projection along the anterior/posterior direction was calculated for each time-point by summing the voxels in the left/right direction across the spinal cord.The anterior/posterior shift was determined by calculating the maximum of the cross-correlation of the projection at each time-point with the first time-point.The shift was then applied to the image data in a block-circular manner.
Only shifts by an integer number of voxels were applied to avoid the need for an extra interpolation step.This step was performed prior to any other pre-processing.
For all functional data, brainstem structures were separated from cervical volumes at the level of the odontoid process.Subsequently, spinal cord functional data were motion-corrected for x-and y-translations using an in-house implementation of AFNI's 3dWarpDrive following the steps in the Neptune Toolbox (https://neptunetoolbox. com/).Motion-corrected data were smoothed with an in-plane 2D Gaussian kernel with full width at half maximum of 2 mm using a custom in-house script relying on tools from AFNI and FSL, and bandpass filtered (0.01-0.1 Hz) using fslmaths (part of FSL).
Warping parameters for spatial normalisation were determined by segmenting and registering the functional data to the PAM50 spinal cord template (De Leener et al., 2018), via an intermediary subjectspecific T2-weighted 3D volume.Specifically, sct_deepseg_sc (Gros et al., 2019) was used to segment the cord from the cerebrospinal fluid (CSF) on motion-corrected functional data and on T2-weighted structural image (sct_propseg (De Leener et al., 2014) was used for one participant's T2-weighted data where sct_deepseg_sc algorithm failed to detect the cord).Manual intervention was needed for accurate segmentation of functional data and was performed in FSLeyes (McCarthy, 2022).Warping parameters for registration of functional data to the PAM50 template were created by combining warp parameters from (1) registering structural T2-weighted image to functional data utilising manually created disc labels on both images and (2) registering the segmented cord from the T2-weighted image to the PAM50 T2-weighted template via sct_register_to_template (De Leener et al., 2018).These warps were applied to pre-processed and filtered functional data via sct_register_multimodal (De Leener et al., 2018) and all subsequently described analyses were carried out in PAM50 template space.Inverse warp parameters obtained from these steps were used to transform PAM50 template CSF and white matter masks to participant functional space which were used in the physiological denoising step described below.
The physiological noise modelling (PNM) toolbox (Brooks et al., 2008) was used to generate 33 slice-specific regressors accounting for physiological noise based on cardiac and respiratory traces, and CSF signal.A bandpass filter (identical to that used on the functional data, 0.01-0.1 Hz) was applied to nuisance regressors (those generated by the PNM and motion regressors obtained from motion correction as described above) to avoid reintroducing noise into the timeseries (Bright et al., 2017).Regression of physiological noise (cardiac and respiratory), CSF and white matter signal, and motion parameters, along with pre-whitening using FILM were performed in FEAT.The smoothed and filtered data (i.e., the residuals from the previous step) were used for subsequent analyses.

| Temporal signal-to-noise ratio
tSNR was calculated on minimally processed resting-state data to avoid artificially inflating the measure.The data had undergone motion correction only (as described above), to remove the timecourse variability associated with in-scan motion and enable creating subject-specific spinal cord masks (see detailed description of steps taken in generating cord masks above).tSNR maps were created by dividing the mean functional image by its standard deviation.Mean tSNR was extracted for all cervical segmental levels (C1-C8) using subject-specific cord masks and for segmental levels C5-C8 using probabilistic segmental masks from the PAM50 atlas (De Leener et al., 2018) warped to subject-space (binarised and thresholded at 30% likelihood of belonging to that spinal level).tSNR was extracted for all complete datasets (complete restingstate acquisition on both MRI sessions, that is, 28 participants/56 resting-state acquisitions) that passed all other quality assurance steps (see Figure 1 for details).Since there are no established guidelines on cut-offs for inclusion based on data quality in spinal fMRI and given that poor tSNR was associated with signal dropout, high level of motion, and geometric distortions, we opted for a minimum tSNR of 20 to ensure reliability estimates were not affected by poor data quality (see Figure S1).Consequently, five participants (i.e., 10 restingstate acquisitions) were excluded from all further analyses due to low mean tSNR across the whole cervical cord (<20) on at least one study session.Spatial extent of resting-state networks at group level was assessed using randomise (Winkler et al., 2014) with threshold-free cluster enhancement (5000 permutations, p < .003[p = .05,Bonferroni corrected for 16 individual seed regions]) and a cord mask obtained from the PAM50 atlas.Permutation-based approach was chosen for this analysis given that spinal cord fMRI data do not fulfil the assumptions necessary for the implementation of random field theory, due to the anisometric shape of the cord.This analysis was performed separately for each session.

| Seed-to-seed connectivity
In addition to the above preregistered seed-to-voxel analysis, a more focused seed-to-seed correlation analysis was performed to assess the strength of connections between regions using Python 3.8.4.
Pearson correlations were computed between each pair of seed regions at subject level using numpy.corrcoeffunction (Harris et al., 2020).The resultant correlation coefficients were Ztransformed using numpy.arctanh(Harris et al., 2020).Statistical significance at group-level was assessed using a one-sample t-test calculated using scipy.stats.ttest_1samp(Virtanen et al., 2020)

| Intraclass correlation coefficient
To systematically evaluate the test-retest performance, inter-session reliability was estimated using: where BMS is the between-target mean squares, EMS is the error mean squares, and k is the number of repeated sessions (as described in Caceres et al., 2009).
ICC values were calculated for each voxel (i.e., voxelwise) using the locally developed ICC toolbox (Caceres et al., 2009)  network was obtained using a one-sample t-test of the first session with a voxelwise t-statistic threshold of 3.5 (equivalent to p = .001)conducted in SPM8 as per the original manuscript describing this approach (Caceres et al., 2009).ICC(3,1) was calculated for each COPE separately.Median ICC values are reported, defined as the reliability measure obtained from the median of the ICC distributions within regions.In addition to this pre-registered approach, additional ICC values were also computed to provide a more detailed understanding of the test-retest reliability of spinal restingstate data.
ICC(3,1) of the mean activation within a network was also computed.Mean signal was extracted from group-level maps obtained from randomise (as described above) using a binarized mask defined from the activation map of Session 1.
Additionally, ICC(3,1) values were calculated on the subject-level Z-scores describing each of the connections in the seed-to-seed analysis.
Finally, ICC(3,1) was calculated for tSNR values extracted from the whole cord and from segmental levels C5-C8 (see below).SPSS v28.0.1.1 with Python3 integration was used to calculate ICC values for mean activation within the network, seed-to-seed connectivities, and tSNR.
Following previous recommendations (Fleiss et al., 2013), ICC values will be categorised accordingly: <0.4 as poor, 0.4-0.59 as fair, DSC ranges from 0 to 1 with higher values indicating better overlap between two sets/maps.A value of 1 would thus correspond to perfect overlap, while a value of 0 would correspond to no overlap.

| Temporal signal-to-noise ratio
To assess signal quality, tSNR was extracted from minimally processed data (motion correction only) for all complete datasets (i.e., prior to excluding participants with mean tSNR across the whole cord <20).
Mean tSNR for the whole cord and segmental levels C5-C6 across sessions are given in Table 1.

| Seed-to-voxel connectivity
To assess the spatial extent of cervical spinal resting-state networks, we estimated seed-to-voxel connectivity maps for each subject and session.This section describes the results of the analysis of data from Session 1 (the corresponding analysis of session 2 data is provided in the Supplementary Material S1).For each seed and segmental level, we observed a statistically significant organisation of spinal restingstate networks ( p < .003).Each seed gave rise to a connectivity pattern that was largely confined to the segment, with sparser betweensegment connections (Figure 3, full spatial maps can be accessed on https://identifiers.org/neurovault.collection:13616).While the spatial extent of clusters was similar across the four quadrants of each segment, we qualitatively observed a dorsal bias in functional connectivity of dorsal seeds and a ventral bias in functional connectivity of ventral seeds.Qualitatively, clusters estimated from Session 2 data had highly similar spatial extent (see Supplementary Material S1 for results of Session 2 data analysis and Figure S2 for overlap between Sessions 1 and 2 maps).

| Seed-to-seed connectivity
To assess the strength of functional connections between horns of the cervical spinal cord, we conducted seed-to-seed correlations T A B L E 1 Temporal signal-to-noise ratio (tSNR) across whole cord and within spinal segmental levels C5-C8 for data acquired on magnetic resonance imaging (MRI) Sessions 1 and 2.   2 and for each of the seed-to-seed connectivities in Figure 5.
On average, voxelwise assessments of ICC in the entire cord Finally, to assess the test-retest reliability of signal quality, ICC values were calculated for tSNR.Across the whole cervical spinal cord captured by our data, tSNR reliability was good (ICC = 0.7).Within segmental levels, tSNR reliability was good for segments C6 (ICC = 0.7), C7 (ICC = 0.6), and C8 (ICC = 0.7), and fair for segment C5 (ICC = 0.5).

| Dice similarity coefficient
DSC assessed the spatial agreement of group-and subject-level resting-state maps between the two sessions.DSC for each network at group and subject level is given in Table 3.
Near-perfect agreement was observed in group-level maps (mean DSC = 0.88 ± 0.03) and good agreement was seen in subject-level maps (mean DSC = 0.67 ± 0.11).

| DISCUSSION
This study investigated cervical spinal cord resting-state networks and their test-retest reliability using a novel acquisition method.In The first aim of this study was to quantify the spatial extent of spinal cervical resting-state networks.Our findings of dorsal and ventral bias in the spatial representations of resting-state networks in the cervical spinal cord are in line with our predictions and complement previous investigations characterising the intrinsic activity of the spinal cord (Barry et al., 2014(Barry et al., , 2016;;Eippert et al., 2017;Kong et al., 2014;Vahdat et al., 2020).In fact, the emergence of distinct sensory (dorsal) and motor (ventral) networks within the cervical 1) for each pair of seed regions.(b) Average ICC between each type of region pair (DH-DH, VH-VH, and DH-VH) grouped by distance between spinal segmental levels (within level: C5-C5, C6-C6, C7-C7, and C8-C8, between neighbouring levels: C5-C6, C6-C7, and C7-C8, and between distant levels: C5-C7, C5-C8, and C6-C8).DH, dorsal horn; ICC, intraclass correlation coefficient; L, left; R, right; VH, ventral horn.
spinal cord has been demonstrated with several different analytical approaches, including data-driven independent component analysis (Kong et al., 2014;San Emeterio Nateras et al., 2016) and hypothesisdriven temporal correlation between ROIs (Barry et al., 2014(Barry et al., , 2016;;Eippert et al., 2017).Further, these networks have been observed both at conventional MR field strength (3 T; Eippert et al., 2017;Kong et al., 2014;Liu et al., 2016;San Emeterio Nateras et al., 2016;Vahdat et al., 2020) and at ultra-high field (7 T; Barry et al., 2014Barry et al., , 2016)).Here, we further confirm the presence of the previously reported dorsodorsal and ventro-ventral cross-talk (Barry et al., 2014(Barry et al., , 2016;;Eippert et al., 2017) with seed-to-seed correlations and further show the emergence of unilateral dorsal and ventral networks (Kong et al., 2014) with seed-to-voxel analyses.Our findings support the notion that these networks reflect intrinsic spinal activity, which mirrors the functional neuroanatomy of the spinal cord.
In addition to the distinct dorsal and ventral networks, we observed a strong within-hemicord (i.e., ipsilateral) connectivity between dorsal and ventral horns of the cervical spinal cord.This is in contrast to previous reports of weak dorsal-ventral connectivity within the hemicord (Barry et al., 2014;Eippert et al., 2017).Nonetheless, strong within-hemicord connectivity between dorsal and ventral horns was observed in non-human primates (Chen et al., 2015) and in one study of a small group of healthy adult volunteers (Weber et al., 2018).Furthermore, dorsal-ventral connectivity was also observed in some participants at ultra-high field, however, these results were not consistent and did not emerge at group level (Barry et al., 2014).Dorsal-ventral connectivity may represent a distinct sensory-motor spinal network, which could support motor reflexes and other more lateralised processing (Chen et al., 2015;Harrison et al., 2021).Indeed, anatomical spinal circuits that connect ipsilateral dorsal and ventral horns, including the monosynaptic stretch reflex and nociceptive withdrawal reflex, are well documented (Pierrot-Deseilligny & Burke, 2012).Nonetheless, given the close proximity of ipsilateral dorsal and ventral horn seeds ($1 mm), relatively large inplane voxel size (1.875 mm Â 1.875 mm), and data processing steps on the detectability of within-hemicord connectivity, further study is needed to establish whether these anatomical circuits contribute to a tertiary spinal resting-state network.Consequently, the current study alone should not be used to support the existence of dorsal-ventral functional connectivity in the human cervical spinal cord.
Similar to previous studies (Kinany et al., 2020;Kong et al., 2014;San Emeterio Nateras et al., 2016), we observed that spinal restingstate networks were largely limited in the rostro-caudal extent, mirroring the segmental organisation of the spinal cord.However, we also observed sparse between-segment connections.Intersegmental connectivity has been reported previously (Eippert et al., 2017; T A B L E 3 Group-level and mean subject-level DSC for each resting-state network. Abbreviations: DH, dorsal horn; DSC, dice similarity coefficient; L, left; R, right; VH, ventral horn.Harita & Stroman, 2017;Ioachim et al., 2019;San Emeterio Nateras et al., 2016;Vahdat et al., 2020) and is thought to reflect ascending sensory and descending motor pathways.In line with our findings, others have reported a decrease of connectivity beyond one vertebral level (Harita & Stroman, 2017;Liu et al., 2016;San Emeterio Nateras et al., 2016;Weber et al., 2018) and, in some cases, weak anti-correlation between regions of different segmental levels (Kinany et al., 2020;Kong et al., 2014).This pattern of results was also observed in this study, with an anti-correlation between right ventral horn at C8 and both ipsilateral and contralateral dorsal horn of segment C6.Such negative relationships may reflect processes related to intersegmental inhibition, perhaps contributing to reflexive actions, proprioception, and nociception (Friesen & Cang, 2001;McBain et al., 2016).
Our second aim was to assess whether cervical spinal restingstate networks could be reliably detected across different scanning sessions.The mixed findings observed in our reliability analysis are in contrast to our predictions and previous reports of good and fair reliability of resting-state connections in the cervical spinal cord, albeit when tested within the same scanning session (Barry et al., 2016;Kaptan et al., 2022;Liu et al., 2016).Test-retest reliability is known to reduce with longer lag between sessions across various contexts (Calamia et al., 2013;Duff, 2012), including brain fMRI (Bennett & Miller, 2010, 2013) and specifically resting-state paradigms (Niu et al., 2020;Yang et al., 2022).Changes related to development, aging, learning, and attention, along with other neuroplastic processes likely underpin the biological reasons for poorer reliability in the long term (Bennett & Miller, 2010, 2013).Furthermore, in cerebral fMRI, the highest reliability is usually achieved in data collected within the same scanning session (Shehzad et al., 2009;Wang et al., 2013), which likely reflects additional impact of scanner characteristics (An et al., 2017).
Given that spinal cord fMRI acquisition is considerably more challenging than brain fMRI, with greater impact of baseline physiology and field inhomogeneities related to surrounding tissues, lower intersession test-retest estimates are to be expected.
In recent years, the reliability and reproducibility of neuroimaging results more broadly has been brought into question (Botvinik-Nezer et al., 2020;Poldrack et al., 2017), with largely mixed evidence of reliability across both tasks (Elliott et al., 2020;Kragel et al., 2021) and resting-state brain fMRI (Noble et al., 2019;Noble, Spann, et al., 2017).In fact, many estimates of brain resting-state connectivity achieve ICC values within the poor range (<0.4) across different resting-state metrics, including voxelwise and region-to-region connectivity (Noble et al., 2019;Noble, Scheinost, et al., 2017).Consequently, the test-retest estimates observed here for spinal cord resting-state networks are similar to those routinely observed in the brain.Furthermore, the spatial extents of these networks were similar across sessions.This suggests that while intensity changes in individual voxels and clusters may differ between sessions, the networks are characterised by a consistent spatial representation over time.
Aside from psychological influences, several factors have been identified, that contribute to low fMRI reliability, including poor tSNR (Bennett & Miller, 2010;Raemaekers et al., 2007), sub-optimal data processing choices (Barry et al., 2016), and confounding effects of motion and/or other non-specific signal changes (Gorgolewski et al., 2013;Noble et al., 2019).The inherent challenges of acquiring spinal cord fMRI recordings, likely result in a compound effect of these factors, which may lead to somewhat lower test-retest reliability estimates than those of brain fMRI (Barry et al., 2016).The continued efforts to improve the quality of spinal cord recordings and finetune pre-processing pipelines will likely help to increase the reliability of spinal fMRI.
Nonetheless, it is important to recognise that high reliability does not always reflect data validity.For instance, it has been observed that correction for artefactual signal, such as motion and physiological noise, can lower test-retest reliability in the brain (Birn et al., 2014;Lipp et al., 2014;Noble et al., 2019;Noble, Spann, et al., 2017) and spinal cord (Kaptan et al., 2022).This likely represents more systematic properties of noise within the data (e.g., regular repetition of cardiac and/or respiratory processes, CSF pulsation leading to cord motion) compared with intrinsic activity within the cord, which may be characterised by more dynamic processes (Kinany et al., 2020).This is further supported by our observation of good reliability of the average tSNR of minimally processed data contrasting with the lower reliability of resting-state networks estimated from the same data.
Consequently, it is vital to consider data reliability and validity together and avoid data processing choices which, while boosting reliability, might have an undue effect on validity.
Most spinal cord fMRI studies use z-shimming alone (Eippert et al., 2017;Kong et al., 2014;Vahdat et al., 2020).While not a primary intention of our study, we did observe that the y-shimming (and to a lesser extent x-shimming) gradients did provide additional signal recovery (Tsivaka et al., 2023).One previous study has also reported dynamic x-, y-, and z-shimming (Islam et al., 2019), which differed from our implementation by applying the linear shimming gradients throughout the EPI acquisition for each slice rather than as gradient lobes.Additionally, we used spectral-spatial excitation pulses for our fMRI acquisition.Since these are designed to only excite water, no additional fat saturation pulses were required, which would have increased the TR needed to acquire images from 38 slices (or reduced the number of slices that could be acquired with the same TR).To date, spinal fMRI has been predominately implemented on Siemens scanners with only few exceptions (e.g., Islam et al., 2019).Our acquisition sequence uses a GE scanner platform and thus provides an alternative to the typically used Siemens-based methods.
The acquisition method described here achieved superior signal quality in comparison to reports describing other sequences used in the field to date, reaching an average tSNR of 26 across scanning sessions.This represents large gains over previously described methods, where average tSNR of spinal EPI data at 3 T typically ranges from 5 to 20 (Barry et al., 2018;Eippert et al., 2017;Kinany, Pirondini, Mattera, et al., 2022;Oliva et al., 2022;Powers et al., 2018).This boost in signal quality may be partly due to the larger in-plane voxel size used in this study (1.875 mm Â 1.875 mm compared with 1 mm Â 1 mm typically used elsewhere; Eippert et al., 2017;Harita & Stroman, 2017;Kong et al., 2014;Liu et al., 2016;San Emeterio Nateras et al., 2016).
Aside from differences in voxel sizes, compared with brain fMRI, the low tSNR of spinal fMRI data is additionally driven by baseline physiology inducing spinal cord motion and CSF pulsation (Piché et al., 2009), and susceptibility artefacts arising from the distinct magnetic susceptibility profiles of surrounding tissues, resulting in signal dropout and image distortions (Saritas et al., 2014).While the tSNR achieved by our acquisition sequence remains lower than that of a typical brain EPI (tSNR of $50 when calculated on minimally processed data; Murphy et al., 2007;Oliva et al., 2022) Barry et al., 2014Barry et al., , 2016) ) and those not including spatial smoothing (Eippert et al., 2017;Kong et al., 2014), including in this study (Supplementary Material S1), suggests that these were unlikely confounds in our data.
It is also important to consider that current best practices for spinal cord fMRI data modelling rely on assumptions that have been validated for cerebral fMRI but not studied in detail in the cord.For instance, early evidence suggests that frequencies higher than the conventional 0.08 Hz cut-off used for brain fMRI (Biswal et al., 1995), may be important drivers of spinal cord signalling (Barry et al., 2016).
Here, we used bandpass filtering of 0.01-0.1 Hz to allow for those higher frequencies, while keeping within the bounds of BOLD-validated frequency distribution.Nevertheless, the neurophysiological mechanisms underpinning assumptions crucial for fMRI data modelling, such as BOLD frequency distribution and haemodynamic response, require further study and validation in the cord.
Although we aimed to obtain 30 complete datasets, and while 37 participants completed one scanning session and 32 completed both sessions, the challenges associated with spinal cord fMRI acquisition and resultant data quality concerns meant that our final sample size was reduced to 23.Longer scanning time due to shimming optimisation, an additional anterior array coil resting on the participant's neck and chest, head and neck positioning minimising neck curvature, and the use of external physiology monitoring equipment likely contributed to the discomfort associated with scanning, increased attrition rate, and led to higher in-scan motion.Further data exclusion was related to low tSNR and signal dropout, some of which may be a result of individual differences in the anatomy of surrounding tissues.High data attrition may be an inevitable attribute of spinal cord fMRI studies and needs to be accounted for during study design and recruitment.
Finally, our study investigated the test-retest reliability of cervical spinal resting-state networks across two separate sessions separated by several days or weeks, while previous studies looked at within-session reliability (Barry et al., 2016;Kaptan et al., 2022;Liu et al., 2016).
The inter-scan interval was variable across participants (1-84 days) due to scanner availability and synchronising the sessions with female participants' menstrual cycles for the purposes of pain-related assessments performed as part of other aims of this study (Kowalczyk et al., 2021).Consequently, while participants adhered to systematic pre-scanning lifestyle guidelines (see Section 2) and the time of the day of the sessions was consistent for each individual to minimised external sources of variability, we cannot ascertain that the variability in the inter-scan interval did not act as a confound.A full characterisation of spinal cord fMRI reliability demands acquiring recordings from the same participants within the same session, as well as over days, weeks, months, and possibly years.Furthermore, combining recordings from the same subject across several sessions has been hypothesised to improve reliability alongside validity (Noble, Spann, et al., 2017).Such efforts in spinal cord fMRI may help to better understand the neurofunctional characteristics of spinal cord restingstate networks.

| CONCLUSIONS
In this study, we demonstrate functional connectivity relationships in dorsal and ventral regions of the cervical cord using a novel, custom acquisition method implemented on a GE platform.Importantly, our findings are in agreement with the known neuroanatomical and neurofunctional organisation of the spinal cord.Although the test-retest reliability of these networks was mixed, their spatial extent was highly reproducible across sessions, suggesting that these networks are characterised by a consistent spatial representation over time.

Conceptualization
and cervical spinal cord to vertebral level T1 (TR = 2.5 s, echo time [TE] = 120 ms, echo train length = 78, flip angle = 90 , field of view [FOV] = 300 mm, acquisition matrix = 320 Â 320, slice thickness = 0.8 mm).This acquisition was based on Cohen-Adad et al. (2021) with the FOV increased to 300 mm.Functional data were acquired over 38 sequential slices in descending order (slice thickness = 4 mm, slice gap = 1 mm), with the inferior-most slices prescribed at vertebral level T1 and covering the whole cervical spinal cord and the brainstem (TR = 2.5 s, TE = 30 ms, flip angle = 90 , ASSET factor = 2, FOV = 180 mm, acquisition matrix = 96 Â 96, reconstruction matrix = 128 Â 128, inplane voxel size = 1.875 mm Â 1.875 mm).Static 0th, 1st, and 2nd order shims were optimised.A spectral-spatial excitation pulse was used to excite only tissue water.Slice-specific linear shims were implemented by adding 0.6 ms duration x-, y-, and z-gradient lobes after the excitation pulse.High-order shimming and x, y, and zshimming were optimised over elliptical regions of interest (ROIs) covering the brain (for slices including the brain) or cord (for slices including the spinal cord).ROIs were drawn manually by the researcher (OSK or SM).Altogether, shim optimisation took approximately 20 min.To maintain consistency and avoid potential systematic differences in ROI drawing affecting test-retest estimates, the same researcher drew ROIs for both MRI sessions within participant.
Assessment of resting-state networks 2.6.1 | Definition of seed regions Seed regions were derived from the PAM50 atlas(De Leener et al., 2018) and corresponded to the four grey matter horns (ventral/ dorsal and left/right) of fifth, sixth, seventh, and eighth segmental levels.We focus on segmental levels C5-C8 given the dermatome projections to the upper limbs and thus their importance in sensorimotor investigations.To obtain these masks we: (1) thresholded the mask of each horn (left/right, dorsal/ventral) at 50% likelihood of belonging to that grey matter horn and binarised it, (2) thresholded probabilistic segmental level (spinal levels C5-C8) masks at 30% to avoid overlap between segments, (3) multiplied each horn mask by each segmental level mask.This resulted in 16 individual masks for seed regions reflecting left/right and dorsal/ventral horns at segmental levels C5, C6, C7, and C8 (Figure2).
2.6.2 | Seed-to-voxel connectivityMean timecourses extracted from these seed regions were used to estimate functional connectivity maps between the seed and each voxel in the cervical cord.For each subject, to assess both within-and between-segment connectivity all four seeds' mean timecourses (left dorsal horn [L DH], right dorsal horn [R DH], left ventral horn [L VH], right ventral horn [R VH]) for a given segmental level (C5, C6, C7, and C8) were included in a single model estimated by FEAT (Z = 3.1, p < .05).Consequently, a total of four models per subject, per session were run.Contrast of parameter estimates (COPE) images from this stage was registered to PAM50 space using warp parameters generated during pre-processing (see above).
. A positive false discovery rate (FDR) was used to account for multiple comparisons (thresholded at p < .05,implemented with statsmodels.stats.multitest.fdrcorrection;Seabold & Perktold, 2010).FDR correction was chosen following convention used in network-based analyses and to avoid excessively conservative thresholding.The analysis presented in the main text of the manuscript used data acquired on Session 1 (see Supplementary Material S1 for corresponding analysis of data acquired on Session 2).
running in Matlab version 9.5.0 (Mathworks Inc.).Reliability was calculated for the whole cord and the complete activation network.The activation F I G U R E 2 Seed regions used in assessments of spinal cord resting-state networks.A total of 16 seeds were derived from the PAM50 atlas, corresponding to the four grey matter horns of the cord at spinal segmental levels C5 (red), C6 (blue), C7 (green), and C8 (yellow).
between each pair of seed regions on data acquired during Session 1 (for results of the same analysis performed on Session 2 data, see Supplementary Material S1).A correlation matrix depicting cervical spinal cord connections is shown in Figure 4. On average, within segment, the strongest statistically significant positive correlations were observed within hemicord (i.e., left DH-VH and right DH-VH), followed by VH-VH and DH-DH connections, and DH-VH connections between hemicords (i.e., left DH-right VH, right DH-left VH).Weaker but statistically significant positive correlations were also observed between neighbouring segments, including DH-DH, VH-VH, as well as within and between hemicords.Finally, negative correlations were observed between the right VH of segment C8 and both left and right DH of segment C6.A similar pattern of results was observed in the analysis of data acquired during Session 2. We further investigated whether these findings were produced by mixing of seed timecourses due to spatial smoothing, repeating our analysis on unsmoothed data.F I G U R E 3 Resting-state networks were obtained from seed-to-voxel connectivity analysis for each of the four horns (ventral/dorsal and left/ right) of segmental levels C5-C8 (data acquired on MRI Session 1).Axial slices are marked with the z MNI coordinate.Each resting-state map was thresholded at p < .003( p = .05,Bonferroni corrected for 16 individual seed regions).F I G U R E 4 Seed-to-seed correlation matrix displaying z-transformed Pearson r.DH, dorsal horn; L, left; VH, ventral horn; R, right.*p < .05,**p < .001.The resultant correlation matrices (detailed in Supplementary Material S1) largely resembled those obtained from the analysis of smoothed data.3.3 | Test-retest reliability 3.3.1 | Intraclass correlation coefficient ICC(3,1) was used to examine the test-retest reliability of cervical resting-state networks.ICC values for each resting-state network derived from seed-to-voxel connectivity analysis are given in Table

(
mean across networks ICC = <0.1 ± <0.1) and within the activation network defined based on MRI session 1 (mean across networks ICC = 0.1 ± <0.1) showed poor reliability across resting-state networks.ICCs for mean activation within each resting-state network showed better but still poor reliability (mean across networks ICC = 0.3 ± 0.2).Nonetheless, more variability in ICC values was observed, with some networks reaching fair (left and right DH networks at level C5 and right VH networks at levels C7 and C8) and good reliability (left VH networks at levels C5 and C6).ICCs for connection strength across pairs of seed regions were variable.ICCs for a large portion of seed pairs (84%) were poor, however some reached fair (14%) and good (2%) levels.Fair and good ICCs were observed for connections both within and between spinal segmental levels and largely reflected either within (i.e., left DH-VH or right DH-VH) or between hemicord connectivity (i.e., left DH-right VH or right DH-left VH).